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Abstract 


We  compare  an  inverse  problem  approach  to  parameter  estimation  with  homogenization  techniques  for 
characterizing  the  electrical  response  of  composite  dielectric  materials  in  the  time  domain.  We  first  consider 
an  homogenization  method,  based  on  the  periodic  unfolding  method,  to  identify  the  dielectric  response  of  a 
complex  material  with  heterogeneous  micro-structures  which  are  described  by  spatially  periodic  parameters. 
We  also  consider  electromagnetic  interrogation  problems  for  complex  materials  assuming  multiple  polariza¬ 
tion  mechanisms  with  distributions  of  parameters.  An  inverse  problem  formulation  is  devised  to  determine 
effective  polarization  parameters  specific  to  the  interrogation  problem.  We  compare  the  results  of  these 
two  approaches  with  the  classical  Maxwell-Garnett  mixing  model  and  a  simplified  model  with  a  weighted 
average  of  parameters.  Numerical  results  are  presented  for  a  specific  example  involving  a  mixture  of  ethanol 
and  water  (modeled  with  multiple  Debye  mechanisms).  A  comparison  between  each  approach  is  made  in 
the  frequency  domain  (e.g.,  Cole-Cole  diagrams),  as  well  as  in  the  time  domain  (e.g.,  plots  of  susceptibility 
kernels) . 
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1  Introduction 


In  this  paper,  we  consider  composite  materials  that  are  mixtures  of  linear,  isotropic  and  temporally 
dispersive  dielectrics.  Composite  materials  are  often  used  in  industry  as  they  have  better  prop¬ 
erties  (e.g.,  electrical,  structural)  as  compared  with  their  constituents.  The  dielectric  properties 
of  mixtures  can  differ  significantly  from  the  individual  materials  from  which  they  are  composed. 
We  compare  three  different  time-domain  approaches  to  characterizing  the  electrical  response  of 
composite  dielectrics.  We  also  consider  the  approach  of  using  the  weighted  average  of  the  model 
parameters  from  each  of  the  constituent  materials. 

First  we  consider  an  homogenization  method  to  identify  the  dielectric  response  of  a  material 
presenting  heterogeneous  micro-structures  which  are  described  by  spatially  periodic  parameters. 
When  such  composite  materials  are  subjected  to  electromagnetic  fields  generated  by  currents  of 
varying  frequencies,  and  the  period  of  the  structure  is  small  compared  to  the  wavelength,  the 
spatially  dependent  coefficients  in  Maxwell’s  equations  oscillate  rapidly.  These  spatially  oscillating 
coefficients  are  difficult  to  treat  numerically  in  simulations.  In  the  process  of  homogenization 
the  rapidly  oscillating  coefficients  are  replaced  with  new,  effective  constant  (in  space)  coefficients. 
The  primary  objective  of  homogenization  is  to  replace  a  system  with  periodically  varying  spatial 
coefficients  by  a  limiting  homogeneous  system  that  facilitates  computation.  The  approach  that 
we  take  here  is  based  on  the  periodic  unfolding  method  (and  its  application  to  homogenization  of 
Maxwell’s  equations)  as  presented  in  [1,  6,  8]. 

In  a  second  method  we  consider  a  parameter  estimation  problem  which,  rather  than  producing  a 
limiting  system  with  effective  coefficients,  chooses  effective  coefficients  for  a  simplified  system  that 
best  matches  the  desired  output  from  a  specific  problem  formulation.  In  particular,  we  assume 
an  interrogation  problem  with  a  given  source  signal  and  known  geometry.  The  data  to  be  fit  in  a 
least  squares  sense  is  a  simulation  of  the  complex  material  using  a  model  involving  a  distribution 
of  polarization  mechanisms  as  described  in  [4] . 

We  compare  the  homogenization  and  parameter  identification  technique  with  results  obtained 
using  the  Maxwell-Garnett  mixing  model  [15]  of  a  two-phase  mixture  wherein  the  inclusions  are 
embedded  in  a  matrix  material  in  the  form  of  spheres  (in  3D)  or  circular  disks  (in  2D)  [14].  We 
also  compare  the  various  techniques  mentioned  above  to  a  single  polarization  model  which  uses 
the  weighted  average  of  the  dielectric  parameters  of  the  two  different  constituents  of  the  composite 
material. 

The  outline  of  this  paper  is  as  follows.  In  Section  2,  we  describe  Maxwell’s  equations  and 
corresponding  constitutive  relations  which  govern  the  evolution  of  the  electromagnetic  field  in  a 
temporally  dispersive  media.  This  is  the  framework  for  the  various  techniques  mentioned  above. 

In  Section  3  we  present  an  outline  of  an  homogenization  approach,  based  on  the  periodic  un¬ 
folding  method,  for  calculating  effective  parameters  of  composite  media.  Using  this  approach  we 
can  calculate  the  high  frequency  permittivity  of  the  composite  dielectric  as  well  as  a  susceptibility 
function  for  discrete  time  values.  In  this  manner  we  obtain  only  time-domain  data.  However,  for 
ease  of  comparison  with  frequency  domain  homogenization  methods,  we  can  calculate  the  complex 
permittivity  of  the  mixture  by  performing  a  discrete  Fourier  transformation  of  the  susceptibility 
function  data. 

In  Section  4  we  introduce  the  parameter  estimation  method  based  on  an  inverse  problem  formu- 
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lation.  Given  data,  and  a  proposed  model  for  simulating  an  approximation  to  the  data,  the  inverse 
problem  formulation  allows  one  to  determine  the  model  parameters  which  best  match  the  data  in 
the  least  squares  sense.  In  the  case  where  the  proposed  model  uses  a  single  polarization  mechanism 
with  scalar  parameters,  the  inverse  problem  method  produces  effective  model  parameters  similar  to 
an  homogenization  method.  However,  the  approach  is  sufficiently  general  to  return  a  distribution  of 
parameters  and/or  mechanisms  (see  [4]).  Furthermore,  this  parameter  estimation  method  is  capa¬ 
ble  of  using  either  simulated  or  real  data  in  the  inverse  problem.  The  method  makes  no  assumption 
on  the  geometry  of  the  mixture,  rather  it  assumes  a  perfectly  homogeneous  composite.  However, 
as  this  approach  does  use  information  about  the  interrogating  fields,  unlike  the  homogenization 
methods,  it  can  be  tailored  to  specific  nondestructive  evaluation  (NDE)  applications  with  narrow 
frequency  bands. 

The  Maxwell-Garnett  mixing  model  presented  in  Section  5  is  also  an  homogenization  type 
method,  but  in  order  to  distinguish  it  from  the  method  described  in  Section  3,  we  will  refer  to 
it  only  as  a  mixing  rule.  It  can  be  used  to  calculate  an  effective  complex  permittivity  of  the 
composite  in  the  frequency  domain  as  well  as  the  susceptibility  function  values  in  the  time  domain. 
In  the  Maxwell-Garnett  model,  it  is  assumed  that  one  of  the  components  of  the  mixture  exists  as 
spherical  (3D)  or  circular  (2D)  inclusions  inside  a  background  material  (ellipsoidal  inclusions  can 
also  be  considered).  The  homogenization  method  assumes  that  the  composite  has  a  periodic  micro¬ 
structure  where  one  of  the  components  in  the  mixture  exists  as  inclusions  of  a  certain  geometry 
inside  a  background  matrix  made  up  of  the  other  component  of  the  mixture.  However,  unlike  the 
Maxwell-Garnett  rule,  the  geometry  of  the  micro-structure  is  not  limited  to  circular  or  spherical 
geometry.  For  a  proper  comparison,  in  this  paper  we  do  assume  a  circular  geometry  for  the 
micro-structure.  For  both  the  homogenization  method  based  on  periodic  unfolding  as  well  as  the 
Maxwell-Garnett  model,  knowledge  of  the  volume  proportions  of  the  components  of  the  mixture  is 
needed  in  order  to  calculate  effective  parameters  of  this  composite  dielectric. 

In  Section  6,  we  present  comparisons  of  these  three  techniques,  along  with  a  simple  model  in¬ 
volving  the  weighted  average  of  the  dielectric  parameters  from  each  constituent  material.  Time 
domain  comparisons  include  plots  of  the  susceptibility  function,  while  the  frequency  domain  com¬ 
parisons  involve  various  plots  of  the  complex  permittivity.  Finally,  we  conclude  with  a  summary  of 
results  and  a  discussion  of  the  advantages  and  disadvantages  of  the  different  techniques  in  Section 
8. 


2  Maxwell’s  Equations  in  a  Dispersive  Medium 

We  employ  Maxwell’s  equations  for  a  linear  and  isotropic  medium  in  a  form  that  includes  terms 
for  the  electric  polarization  given  by 

<9D 

—  -  VxH  =  J  in  (0,  T)  x  H, 

<9B 

-7—  +  VxE  =  0  in  (0, T)  x  H, 

E  x  n  =  0  on  (0,  T)  x  <9H, 

E(0,x)  =  0,  H(0,x)  =  0  in  H, 


(1) 
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along  with  zero  Gauss  divergence  laws  (V •  D  =  0,  VB  =  0  in  (0,  T)  x  Q).  Here  Q  C  Mn,  n  =  1,2, 3. 
The  vector  valued  functions  E  and  H  represent  the  strengths  of  the  electric  and  magnetic  fields, 
respectively,  while  D  and  B  are  the  electric  and  magnetic  flux  densities,  respectively.  The  source 
current  density  is  given  by  J.  We  assume  that  there  are  no  free  electric  charges  unaccounted  for 
in  the  electric  polarization.  We  use  perfectly  conducting  boundary  conditions  on  the  boundary 
dD,  with  unit  outward  normal  n  and  also  assign  zero  initial  conditions  for  all  the  unknown  fields. 
System  (1)  is  completed  by  constitutive  laws  that  embody  the  behavior  of  the  material  in  response 
to  the  electromagnetic  fields.  These  are  given  in  (0,  T)  x  O  in  the  form 


D(t,  x)  =  e0er(x)E(t,  x)  +  P(t,  x), 
„  B(f,x)  =  /i0H(t,x), 


(2) 


where  eo,  and  /M)  are  the  permittivity  and  the  permeability  of  free  space,  respectively,  er  is  the 
relative  permittivity  of  the  medium  under  investigation  and  P  is  the  media’s  macroscopic  electric 
polarization.  For  the  media  that  is  of  interest  to  us,  we  have  neglected  magnetic  effects  and  we  have 
assumed  zero  conductivity.  To  describe  the  behavior  of  the  media’s  macroscopic  electric  polarization 
P,  we  employ  a  general  integral  representation  model  in  which  the  polarization  explicitly  depends 
on  the  past  history  of  the  electric  field.  This  convolution  model  is  sufficiently  general  to  include 
microscopic  polarization  mechanisms  such  as  dipole  or  orientational  polarization  as  well  as  ionic 
and  electronic  polarization  and  other  frequency  dependent  polarization  mechanisms.  The  resulting 
constitutive  law  can  be  given  in  terms  of  a  susceptibility  kernel  g,  also  known  as  the  dielectric 
response  function  (DRF),  by 

P(i,  x)  =  g  *  E(i,  x)  =  f  g(t  —  s,  x)E(s,  x)  ds.  (3) 

Jo 


2.1  The  Debye  Model  of  Orientational  Polarization 


The  constitutive  law  (2)-(3)  is  sufficiently  general  to  include  models  based  on  differential  equa¬ 
tions  and  systems  of  differential  equations  or  delay  differential  equations  whose  solutions  can  be 
expressed  through  fundamental  solutions  (in  general  variation-of-parameters  representations)-see 
[3]  for  details.  A  number  of  known  polarization  laws  can  be  readily  treated.  The  choice  of  the 
kernel  function 


e0(e8(t,x)  -  foo(t1x))  -t/T(t,x) 
r(t,x) 


(t,  x)  e  (0,  T)  x  H 


(4) 


in  the  dielectric  corresponds  to  the  differential  equation  of  the  Debye  Model  for  orientational  or 
dipolar  polarization  given  by 


r(t,x)Pt(f,x)  +  P(t,x)  =  e0(es(t,x)  -  e^t,  x))E(t,  x), 
D(t,x)  =  e0eoo(t,x)E(t,x)  +  P(t,x). 


(5) 

(6) 


Here  es  is  the  static  relative  permittivity.  The  presence  of  instantaneous  polarization  is  accounted 
for  in  this  case  by  the  coefficient  er  =  in  the  electric  flux  equation.  The  remainder  of  the 
electric  polarization  is  seen  to  be  a  decaying  exponential  with  relaxation  parameter  r,  driven  by 
the  electric  field,  less  the  part  included  in  the  instantaneous  polarization.  This  model  was  first 
proposed  by  Debye  [10]  to  model  the  behavior  of  materials  that  possess  permanent  dipole  moments. 
The  magnitude  of  the  polarization  term  P  represents  the  degree  of  alignment  of  these  individual 
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moments.  The  choice  of  coefficients  in  (5)  gives  a  physical  interpretation  to  es  and  as  the  relative 
permittivities  of  the  medium  in  the  limit  of  the  static  field  and  very  high  frequencies,  respectively. 
In  the  static  case,  we  have  Pj  =  0,  so  that  P  =  eo(e<j  —  Coo)E  and  D  =  eseoE.  For  very  high 
frequencies,  r Pf  dominates  P  so  that  P  «  0  and  D  =  c^coE. 


3  Homogenization  Model  in  Two  Dimensions 


We  assume  that  our  problem  possesses  uniformity  in  the  spatial  direction  x2.  so  that  all  derivatives 
with  respect  to  x2  are  assumed  to  be  zero.  In  this  case  Maxwell’s  equations  decouple  into  two 
different  modes,  the  transverse  electric  (TE)  and  transverse  magnetic  (TM)  modes.  Here,  we 
are  interested  in  the  TE  mode  for  which  the  electric  field  has  two  nontrivial  components,  E  = 
(EXl,  EX3)t,  and  one  nontrivial  component  H  =  HX2  of  the  magnetic  field.  Similarly,  the  flux 
densities  are  D  =  (DX1,  DX3)T,  and  B  =  BX2.  The  system  then  reduces  to  vectors  depending  on 
two  variables,  x\  and  x%. 


Let  H  now  be  a  subset  of  M2  with  x  =  (aq,  X3)  £  M2.  Let  us  define  the  vector  of  fields 


and  the  operator 


Lu(f,x)  = 


£  W1A{0,T;H\n-,R3)), 

(7) 

D (Lx)  \ 

H(t,x)  / 

(8) 

which  from  (2)-(3)  can  be  written  as 


LiEu(f,x)  =  Ai  (x)u(f,x)+  /  CiE(t  —  s,  x)u(s,  x)  ds. 


(9) 


The  two  3x3  coefficient  matrices  in  (9)  are  defined  as 

Ate  = 

with 


"A™ 

O' 

;  CTE  = 

rnTE 

^11 

O' 

_  o2 

Mo. 

_  o2 

0 

A™  (x)  = 


e0er(x)  0 
0  e0er(x) 


CnE(t,x)  = 


g(t,x)  0 
0  g(t,  x) 


(10) 


(11) 


where  in  the  above  definitions  O2  is  the  2x2  matrix  of  zeros.  Next,  we  define  the  Maxwell  curl 
operator  M1E  as 

MTEu(t,x)  =  MTEf  E(*’X)  V  (12) 


H(t,  x) 


where  the  operator  MrE  can  be  expressed  as  the  matrix  operator 


/  0  0  -dX3  \ 

0  0  dxi  \ 

\  dX3  -dXx  0  / 


(13) 


6 


Thus  Maxwell’s  equation  can  be  rewritten  in  the  form 


(i) 

-^-Lteu  =  Mteu  +  J  in  (0,  T)  x  17, 

(ii) 

u(0,  x)  =  0  in  17, 

(14) 

(hi) 

ui(f,x)  x  n(x)  =  0  on  (0,  T)  x  517, 

where  LTE  is  the  operator  associated  with  the  constitutive  law  (9),  and  M1E  is  the  Maxwell 
operator  (12). 

We  assume  that  the  structure  that  occupies  the  domain  17  entails  periodic  micro-structures 
leading  to  matrices  ArE,B1E  and  CIE  with  spatially  oscillatory  coefficients.  Specifically,  we  will 
assume  that  er  and  g  are  rapidly  oscillating  spatial  functions. 


Y“  Cell 


a  | 


Figure  1:  Periodic  composite  material  presenting  a  circular  micro-structure  with  periodicity  a.  The 
figure  shows  a  decreasing  from  left  to  right. 


3.1  The  Homogenized  Solution 

The  outline  presented  in  this  section  is  based  on  results  from  [1,  6].  We  denote  by  Y°  the  reference 
cell  of  the  periodic  structure  that  occupies  17  C  R2  (see  Figure  1).  Let  xjgK2  with  x  =  (aq,  £3) 
and  y  =  (2/1, 2/3).  We  will  use  x  £  R2  for  points  on  the  macro  scale,  and  y  G  M2  for  points  on  the 
micro  scale  (reference  cell).  Since  we  assume  uniformity  in  the  aq-direction,  we  may  have  as  before 
all  derivatives  with  respect  to  a:  2  (or  222  in  the  micro-scale)  zero.  In  this  case  the  homogenized 

solution  u  =  (EX1,  EX3,  HX2)T  for  the  TE  mode  appears  in  a  formal  asymptotic  expansion,  or 

two-scale  expansion  (see  [1]),  given  using  a  corrector  term  u(f,x,y)  by 

Ex!  =  Ex 1  +  dyiu(t,  x,  y)  +  . . . ,  (15a) 

Ex3  =  ex3  +  dy3u(t,  x,  y)  +  . . . ,  (15b) 

H%2  =  HX2,  (15c) 

where  the  electromagnetic  field  u"  =  ( E“] ,  E“3 ,  i7“2 ) T  is  the  solution  to  an  evolution  problem  of 
type  (14)  on  a  domain  17  with  a  micro-structure  of  periodicity  a.  Hence  the  homogenized  electric 


field  solution  for  the  TE  mode  appears  in  the  formal  expansion 

E"  =  E  +  Vy«(f,  X,  y )  +  . . . 
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(16) 

Since  we  are  assuming  uniformity  in  the  X2  direction,  there  is  no  correction  of  the  magnetic  field 
solution. 


The  solution  of  the  first  corrector  term  u  £  VE2,1(0,  T;  H^er{Y]  M)),  from  the  two  scale  ex¬ 
pansion,  depends  on  the  calculation  of  two  corrector  sub-terms  £  i7per(y;M)  and  wk  £ 
VE1,1(0,  T;  HxpjY ;  M)),  which  are  solutions  to  the  corrector  equations 


(i)  J  A™(y)Vy^(y)  •  Vyu(y)  dy  =  -  J '  A^(y)ek  •  Vyu(y)  dy, 

(ii)  J ^  A?’1E(y)  Vy  wk(t,  y)  •  Vyu(y  )dy 

+  [  f  s,y)Vywk{s,y)  ds-X7yv(y)dy 

JY  JO 

=  ~  f  Cn(f,y)  {ek  +  Vyu)^(y)}  ■  Vyv(y)  dy, 

JY 


(17) 


for  all  v  £  i7pCr(y),  k  =  1,2.  The  vectors  ek,  k  =  1,2  are  the  basis  vectors  in  M2  with  ei  =  (1,0)T, 
e-2  =  (0, 1)T .  The  gradient  operator  in  this  case  is  Vy  =  (dyi,  dy3)T .  The  first  corrector  term  u  can 
now  be  calculated  as 

u(t,x.,y)  =  wA(y)E(t,x)  +  (  w(t  -  s,y)E(s,x)  ds,  (18) 

Jo 

where  £  Mlx2  with  components  {te^}2.=1.  Similarly  w  £  Rlx2  with  components  {wk}‘i=  \  ■ 


Once  we  have  solved  for  these  corrector  terms,  we  can  then  construct  the  homogenized  matrices 
from 


(i)  (AnE)fc  =  A^E(y)  {ek  +  Vyu>£( y)}  dy, 

(ii)  {Cjf)k{t)  =  J r  C£E(i,y)  {efc  +  Vyu>^(y)}  dy  +  J ^  Ajf(y)Vyu>fc(i,y)  dy  (19) 

+  [  [  — s,y)Vy«7fc(s,y)  ds  dy, 

Jy  Jo 

where,  (A™)*;  and  (Cj^)k  are  the  fcth  columns  of  the  2x2  matrices  A™  and  C™,  respectively, 
and  the  homogenized  matrices  are  given  by 


ate  = 

1 

H 

O 

,  CTE  = 

"C? iE  o' 

o 

o 

o 

o 

The  corresponding  system  of  equations  in  the  TE  mode  are 

(i)  -^£teu  =  M1Eu  +  jJE  in  (0,  T)  x  12, 

(ii)  u(0,  x)  =  0  in  12, 

(iii)  ui(t,x)  x  n(x)  =  0  on  (0 ,  T)  x  <912, 


(20) 


(21) 


where  u  =  (EX1,  EX3,  HX2)T,  n  =  (nxi,nx 3)T  is  the  unit  outward  normal  vector  to  dfl,  and  the 
operator  £TE  is  defined  as 

£TEu(f ,  x)  =  ATEu(f ,  x)  +  f  CTE(t  —  s)u(s,  x)  ds.  (22) 

Jo 

If  the  initial  conditions  are  nonzero,  additional  corrector  terms  and  corrector  equations  are  involved, 
and  then  there  is  a  supplementary  source  term  J'0  that  should  be  introduced  in  the  right  side  of 
(21, i).  See  [6]  for  more  details. 


3.2  Numerical  Discretization  for  the  Cell  Problem 

The  spatial  discretization  of  the  cell  problem  can  be  done  using  piecewise  bilinear  finite  elements. 
For  the  time  discretization  we  use  an  approach  that  involves  calculating  the  convolution  recursively. 
Since  the  susceptibility  kernel  </(i,x)  is  exponential  in  nature  for  many  materials  of  interest,  we 
can  use  recursion  to  compute  the  discretized  time  convolution  of  the  susceptibility  kernel  with  the 
gradient  of  the  corrector  terms  in  the  corrector  subproblem  (17, ii)  and  the  construction  (19, ii).  The 
discrete  gradient  of  the  corrector  term  is  assumed  to  be  constant  on  each  equally  spaced  interval 
of  time  At.  The  details  of  the  method  summarized  here  may  be  found  in  [1].  A  similar  approach 
known  as  the  recursive  convolution  (RC)  method  has  been  used  to  compute  the  discrete  convolution 
terms  that  appear  in  Maxwell’s  equations  [12,  13].  We  show  in  [1]  that  the  discrete  convolution  of 
all  previous  field  values  of  the  gradient  of  the  corrector  terms  and  the  discrete  susceptibility  function 
can  be  reduced  to  recursive  updating  of  a  single  vector  on  each  element  in  the  finite  element  mesh; 
this  involves  a  matrix  vector  multiplication  at  each  time  step. 

For  the  numerical  experiments  in  this  paper  we  obtain  the  effective  high  frequency  permittivity 
as  well  as  time  plots  of  the  effective  susceptibility  kernels  for  various  dielectric  mixtures.  We  do 
so  by  solving  specific  cell  problems  in  the  reference  cell  Y  =  [0,1]  x  [0,1].  The  simulations  are 
performed  on  a  51  x  51  nodes  finite  element  grid.  In  the  numerical  solution  for  the  homogenized 
model  we  have  approximated  circular  micro-structures  in  a  staircase  fashion.  For  details  see  [1]. 
The  results  of  several  experiments  will  be  presented  in  Section  6. 


4  Parameter  Estimation 


We  next  formulate  an  inverse  problem  for  the  determination  of  parameters  for  a  model  which  de¬ 
scribes  the  dielectric  response  of  complex  materials  in  the  context  of  electromagnetic  interrogation. 
It  is  well  known  that,  in  general,  dielectrics  are  most  accurately  modeled  by  a  distribution  of  pa¬ 
rameters  [7],  although  using  a  distribution  of  parameters  in  practice  is  sometimes  computationally 
intensive.  Thus  there  is  a  need  to  determine  single,  effective  parameters.  The  comparison  of  several 
ways  to  do  this  is  the  precisely  the  subject  of  this  paper. 

In  [4]  it  was  shown  that  the  inverse  problem  for  the  actual  distribution  is  well  posed.  In  fact, 
distributions  of  mechanisms  (e.g.,  mixtures  of  Debye  and  Lorentz  materials)  could  also  be  handled 
in  the  same  framework.  Using  this  theory  as  a  foundation,  we  consider  the  inverse  problem  wherein 
parameters  for  a  simplified  model  are  estimated.  The  most  simplified  case  of  a  single  mechanism 
involving  a  delta  distribution  of  the  parameter  set  results  in  the  determination  of  effective  model  pa¬ 
rameters  as  produced  in  homogenization  techniques.  However,  unlike  homogenization  methods,  the 
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inverse  problem  formulation  has  the  flexibility  to  use  either  real  data  or  simulated  data  obtained  by 
using  a  suitable  forward  solve.  Further,  as  the  electromagnetic  interrogation  of  an  object  generally 
involves  a  relatively  narrow  range  of  frequencies,  the  inverse  problem  based  parameter  estimation 
is  designed  to  essentially  match  the  dielectric  response  of  the  material  in  the  frequency  band  of 
interest.  This  is  unlike  homogenization  methods  which  do  not  take  the  interrogating  frequency  into 
account. 


4.1  Distributions  of  Mechanisms  and  Parameters 

As  before,  we  consider  the  electric  field  E  and  the  magnetic  field  H  to  be  governed  by  Maxwell’s 
equations  in  a  domain  V  =  Do  U  fl  where  P  contains  the  dielectric  and  the  ambient  Slo  is  treated 
as  a  vacuum.  (See  Figure  2  for  a  schematic  of  a  sample  domain.) 


Figure  2:  Snapshot  of  a  simulation  of  a  truncated  sine  wave  partially  reflecting  from  and  partially 
penetrating  a  sample  material.  The  z-axis  represents  a  sample  ID  domain. 


Our  main  focus  in  this  section  is  the  dielectric  polarization  P  which  we  assume  has  the  general 
convolution  form  given  in  (3).  In  every  practical  example  (Debye,  Lorentz,  etc.)  DRF’s  are  para¬ 
meter  dependent  as  well  as  time  (and  possibly  space)  dependent;  we  represent  this  as  g  =  g(t,  x;  u), 
where  typically  v  contains  parameters  such  as  the  high  frequency  limit  dielectric  permittivity  e^, 
the  static  permittivity  es,  and  relaxation  time  r.  An  example  of  an  often-used  DRF  is  the  De¬ 
bye  model  (described  in  Section  2.1).  Allowing  for  a  distribution  F  of  parameters  v  over  some 
admissible  set  Af,  we  generalize  the  polarization  law  (3)  to 

P(t,  x;  F)  =  g  *  E(t,  x)  =  f  f  g(t  —  s,  x;  z/)E(s,  x)dF(u)ds.  (23) 

Jo  J M 

We  expect  to  chose  F  from  (or  from  a  subset  of)  the  space  F  =  *P(jV)  of  all  probability  measures  F 
on  A f.  We  could  further  generalize  (23)  to  allow  for  dielectric  materials  with  multiple  mechanisms 
or  multiple  DRF’s  (i.e.,  heterogeneous  molecular  structures)  by  considering  a  family  of  possible 
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DRF’s  and  distributions  over  this  family  as  described  in  [4] .  For  now  we  will  restrict  our  efforts  to 
single  mechanisms. 


4.2  Inverse  Problem  Formulation 


In  [4]  the  focus  was  to  determine  the  parameters  of  distributions  in  an  inverse  problem  formulation 
in  order  to  accurately  describe  the  electrical  response  of  complex  materials.  In  this  effort  we 
attempt  to  find  a  single  set  of  effective  parameters  that  “best”  describe  this  response  for  a  specific 
interrogating  source.  In  other  words,  we  are  attempting  to  find  the  distribution  F  in  the  subset 
of  T  =  fP(A0  which  corresponds  to  all  single  delta  distributions,  even  though  the  data  may  (and 
does!)  come  from  simulations  involving  more  than  just  a  single  parameter  set. 


Formally,  we  consider  the  Maxwell  system  (1)— (2)  with  polarization  P 
(3).  Let 


z (Lx;  F) 


( E(t, x;  F)\ 
\U (t,x;F)J  ’ 


P(t,  x;  F)  given  by 


with  (t,  x)  — ►  z(t,  x;  F)  mapping  from  (0,  T)  x  H  to  i?6.  We  assume  we  are  given  data  d  = 
{d  i}?=1  corresponding  to  observations  of  Caz(L,  • ;  F).  Here  Ca  denotes  evaluation  of  one  or  more 
components  of  E  or  H  at  an  antenna  {x^}.  We  use  this  data  to  estimate  the  parameter  distribution 
F  in  an  ordinary  least  squares  (OLS)  formulation,  seeking  to  minimize 


J(F)  =  J2MU,-;F)-di\2,  (24) 

i= 1 


over  some  subspace  the  space  F  =  *P(JV)  of  all  probability  measures  F  on  M . 


For  this  particular  effort,  we  are  concerned  with  comparing  the  inverse  problem  based  parameter 
estimation  method  to  other  methods  of  determining  effective  parameters.  Thus,  while  the  synthetic 
data  is  generated  using  the  general  formulation  of  the  polarization  law  (23),  the  forward  solves  in 
the  inverse  problem  will  be  limited  to  a  single  polarization  mechanism  with  a  single  parameter  set 
v,  namely,  equation  (3). 


4.3  Numerical  Results 

As  our  inverse  problem  approach  does  not  depend  on  the  inclusion  geometry,  we  may  consider  a 
1-D  example  as  explained  in  detail  in  [5].  Restricting  to  one  dimension,  and  using  D  =  eE  -\-  P,  we 
can  write  Maxwell’s  equations  in  second  order  form  as 

Mo  eE  +  Mo  InP  ~  E"  =  ~m o<A  in  H  U  =  [0,  b],  (25) 

where  E  is  the  k  or  z  component  of  the  electric  field,  and  e  =  eo(l  +  In(e oo  —  1)).  Here,  Iq 
denotes  the  indicator  function  on  the  dielectric  medium  Q.  The  (now  scalar)  polarization  is  still 
determined  by  equation  (3),  however  with  the  vector  valued  electric  field  E  replaced  by  the  scalar 
E.  We  assume  zero  initial  conditions  for  E  and  P, 


P( 0,  z)  =  E( 0,  z)  =  E( 0,  z)  =  0. 
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Lastly,  we  have  absorbing  boundary  conditions  at  2  =  0  and  reflecting  conditions  at  z  =  b,  namely 

\e-cE'} 


=  0 

J  z= 0 

E(t,b)  =  0, 


(26) 


where  c2  = 

eo«> 

The  numerical  method  employed  involves  a  linear  Finite  Element  Method  for  spatial  discretiza¬ 
tion,  and  a  centered  difference  for  the  temporal  discretization.  For  details,  see  [2]. 

We  choose  to  consider  a  mixture  of  two  Debye  materials,  for  which  the  parameter  identification 
procedure  will  attempt  to  find  parameters  for  a  single  Debye  model  such  that  the  resulting  simu¬ 
lations  best  match,  in  a  least  squares  sense,  the  data  generated  from  the  mixture.  The  substance 
that  is  being  interrogated  is  a  mixture  of  water  and  ethanol  with  volume  fractions  of  80%  —  20%, 
50%  —  50%,  and  20%  —  80%.  For  water,  we  used  the  following  parameters  in  a  Debye  polarization 
model:  T\  =  1.01  x  10-11,  eS;i  =  80.1,  and  1  =  4.9.  The  Debye  parameters  for  ethanol  are 
t-2  =  1.2  x  1CT10,  eS)2  =  25.1,  and  €00,2  =  4.4. 

We  assume  the  interrogating  signal  to  be  a  single  cycle  truncated  sine  curve  with  carrier  fre¬ 
quency  109  Hz.  Note  that  this  frequency  is  below  the  weighted  average  relaxation  frequency 
/  =  1/f,  with  f  =  ctTi  +  (1  —  a) 72  (here  a  is  the  volume  fraction  of  water),  and  therefore  we 
expect  that  the  optimal  single  relaxation  time  should  be  close  to  this  average  value  [4] .  Further,  for 
simplicity  of  computations,  we  take  the  value  of  in  the  single  Debye  model  to  be  the  weighted 
average  of  mixture,  namely  Coo,e//  =  aeoo.i  +  (1  —  aOeoo,2-  This  results  in  a  two  parameter  inverse 
problem  for  r  and  es.  The  interrogating  source  and  the  signal  receiver  are  both  located  at  2  =  0 
(see  Figure  2).  For  this  test  problem,  the  cost  functional  is 

n 

J{r,es)  =  Y \E(U,0;t,c8)  -di\2.  (27) 

i=l 

A  surface  plot  of  each  of  the  cost  functionals  J  is  displayed  in  Figures  3-5  for  the  volume 
fractions  of  80%  —  20%,  50%  —  50%,  and  20%  —  80%,  respectively.  The  left  plot  of  each  figure 
depicts  the  range  of  values  between  the  parameters  for  ethanol  and  the  parameters  for  water.  The 
right  plot  of  each  figure  shows  a  closer  view  of  the  objective  function  near  the  optimal  solution. 
The  initial  starting  point  used  in  the  optimization  method  is  marked  with  an  “O”.  This  value  is 
simply  the  weighted  average  of  r  and  es  over  the  material.  The  resulting  optimum  value  from  the 
optimization  method  is  marked  with  an  “X” . 

Note  that  in  the  80%  —  20%  case  there  is  a  large  positive  peak  centered  at  es  =  60.  We  do  not 
have  to  resort  to  taking  the  absolute  value  of  the  electric  field  to  avoid  this  obstacle  to  minimization, 
as  described  in  [4],  since  our  starting  point  for  the  optimization  method  is  closer  to  the  minimum 
than  this  peak.  For  the  20%— 80%  case  there  is  a  local  minimum  near  es  =  53.  Fortunately,  for  these 
particular  parameters,  the  initial  solution  is  close  enough  to  the  global  minimum  to  not  converge 
to  the  local  minimum.  It  is  important  to  note  however,  that  for  other  choices  of  parameters, 
convergence  to  the  incorrect  minimum  is  a  very  real  concern. 

For  the  optimization  method,  we  use  a  Levenberg-Marquardt  method  with  implicit  filtering.  As 
the  cost  functional  is  significantly  more  sensitive  to  the  es  direction  than  the  r  direction,  taking  a 
Newton  step  in  both  directions  at  the  same  time  generally  results  in  an  increased  value  of  the  cost 
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Figure  3:  The  left  plot  depicts  the  objective  function  for  the  parameter  identification  problem 
versus  the  log  of  r  and  es  using  a  volume  fraction  of  80%  —  20%  water  to  ethanol.  The  right  plot 
presents  a  closer  view  of  the  objective  function,  with  the  initial  value  marked  with  an  “O”  and  the 
optimal  value  marked  with  an  “X” . 


function.  Rather  than  increasing  the  Levenberg-Marquardt  parameter  immediately  (i.e. ,  performing 
a  pull-back  on  the  line-search  which  results  in  a  shorter  step),  we  test  the  possibility  that  taking  just 
one  of  the  directions  will  give  an  actual  reduction  in  the  objective  function.  The  addition  of  these 
two  function  evaluations,  when  the  cost  is  not  decreased  by  the  original  step,  saved  approximately 
70%  of  the  pull-backs  on  line-searches  versus  the  traditional  Levenberg-Marquardt  method,  and 
led  to  faster  convergence.  Table  1  contains  the  initial  values  of  log(r),  es  and  J,  as  well  as  the 
optimized  values  for  each  case  of  differing  volume  fractions  (where  a  is  the  percentage  of  water). 


Table  1:  Initial  values  and  values  resulting  from  Levenberg-Marquardt  optimization  of  the  water- 
ethanol  mixture  using  the  volume  fractions  of  80%  —  20%,  50%  —  50%,  and  20%  —  80%  interrogated 
at  /  =  109  Hz. 


a  =  .8 

log(r) 

es 

J 

Initial 

-10.7807 

69.1 

0.350694 

Optimal 

-10.7933 

67.5529 

0.0122314 

a  =  .5 

log (t) 

es 

J 

Initial 

-10.4582 

52.6 

2.48618 

Optimal 

-10.5759 

48.2837 

0.182077 

a  =  .2 

log(r) 

es 

J 

Initial 

-10.1358 

36.1 

2.92976 

Optimal 

-10.418 

30.1762 

1.54837 

To  see  the  improvement  in  the  fit-to-data  of  the  optimal  value  response  over  that  of  the  weighted 
average  in  matching  the  multiple  Debye  signal,  we  plot  in  Figures  6-8  the  resulting  simulations  and 
the  original  data  computed  using  a  discrete  distribution  of  two  Debye  models,  for  each  of  the  three 
volume  fractions.  The  inset  shows  a  closer  view  of  the  full  signal  reflected  from  the  back  boundary 
(the  only  useful  information  in  determining  the  electrical  response  in  an  experiment). 
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Figure  4:  The  left  plot  depicts  the  objective  function  for  the  parameter  identification  problem 
versus  the  log  of  r  and  es  using  a  volume  fraction  of  50%  —  50%  water  to  ethanol.  The  right  plot 
presents  a  closer  view  of  the  objective  function,  with  the  initial  value  marked  with  an  “O”  and  the 
optimal  value  marked  with  an  “X” . 


It  is  clear  from  each  of  the  figures  that  the  weighted  average  parameters  do  not  result  in  signals 
that  match  the  double  Debye  data.  For  the  cases  80%  —  20%  and  50% — 50%,  the  optimal  parameters 
produce  a  signal  that  is  very  similar  to  the  data,  although  still  distinguishable.  However,  even  for 
these  simple  examples  it  is  apparent  that  no  single  Debye  model  can  produce  a  simulation  exactly 
matching  the  double  Debye  model.  In  fact,  this  conclusion  is  even  more  obvious  in  the  20%  —  80% 
case.  The  optimization  method  returns  this  solution  as  optimum,  and  the  parameter  values  can  be 
visually  verified  by  referring  to  the  surface  plot  5,  but  clearly  the  single  Debye  simulation  does  not 
match  the  multiple  Debye  in  amplitude  nor  phase.  Still,  when  an  approximate  result  is  sufficient 
and  the  savings  in  computational  time  important,  then  the  optimal  parameters  from  an  inverse 
problem  formulation  provide  an  alternative  to  modeling  a  complex  material  using  distributions  of 
parameters  and  mechanisms. 


5  The  Maxwell  Garnett  Model 


We  will  compare  the  homogenization  method  and  the  parameter  identification  technique  presented 
above  to  the  Maxwell-Garnett  (MG)  model  of  a  two-phase  mixture  wherein  the  inclusions  are 
embedded  in  the  environment  in  the  form  of  spheres  (in  3D)  or  circular  disks  (in  2D).  As  in  the 
case  of  the  other  techniques,  we  assume  that  both  the  temporally  dispersive  materials  are  isotropic, 
homogeneous  and  non-magnetic. 

A  direct  time-domain  approach  that  has  been  developed  in  [14]  expresses  the  effective  permit¬ 
tivity  of  a  mixture  as  an  operator  form  MG  expression  that  includes  the  susceptibility  kernels  and 
high  frequency  responses  of  the  component  materials.  The  evaluation  of  the  effective  permittivity 
operator  requires  the  calculation  of  convolutions  and  operator  inverses.  The  inverse  of  the  convolu¬ 
tion  operators  encountered  in  these  problems  can  be  solved  from  Volterra  equations  of  the  second 
kind,  which  have  unique  and  well  behaved  solutions.  The  advantage  of  the  time-domain  formulation 
is  that  the  results  are  more  intuitive  physically  than  those  in  the  frequency  domain.  When  both 
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Figure  5:  The  left  plot  depicts  the  objective  function  for  the  parameter  identification  problem 
versus  the  log  of  r  and  es  using  a  volume  fraction  of  20%  —  80%  water  to  ethanol.  The  right  plot 
presents  a  closer  view  of  the  objective  function,  with  the  initial  value  marked  with  an  “O”  and  the 
optimal  value  marked  with  an  “X” . 


Figure  6:  Comparison  of  multiple  Debye  data  to  the  forward  simulation  using  the  weighted  average 
Debye  parameter  values  and  the  forward  simulation  using  optimal  values  from  the  single  Debye 
inverse  problem  using  a  volume  fraction  of  80%  —  20%. 


constituent  materials  are  dispersive,  the  dispersion  of  the  mixture  is  in  general  more  complicated 
than  that  of  the  inclusion.  Even  for  the  simplest  case  of  Debye  inclusions  in  a  Debye  background, 
the  effective  complex  permittivity  is  neither  of  Debye  nor  of  Lorentz  form,  but  contains  higher 
powers  of  the  frequency  to  in  both  the  numerator  and  the  denominator. 

For  the  examples  of  an  80%  water  and  20%  alcohol  mixture,  and  a  20%  —  80%  mixture,  explicit 
expressions  for  the  effective  susceptibility  kernel  are  given  in  Sec  5.B.3  of  [14].  The  analysis  in 
[14]  is  restricted  to  the  case  of  spherical  (circular  in  2D)  inclusions.  The  homogenization  technique 
presented  in  Section  3,  however,  does  not  have  this  restriction. 


15 


Figure  7:  Comparison  of  multiple  Debye  data  to  the  forward  simulation  using  the  weighted  average 
Debye  parameter  values  and  the  forward  simulation  using  optimal  values  from  the  single  Debye 
inverse  problem  using  a  volume  fraction  of  50%  —  50%. 
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Figure  8:  Comparison  of  multiple  Debye  data  to  the  forward  simulation  using  the  weighted  average 
Debye  parameter  values  and  the  forward  simulation  using  optimal  values  from  the  single  Debye 
inverse  problem  using  a  volume  fraction  of  20%  —  80%. 


6  Comparison  of  Methods 


We  compare  the  time  domain  and  frequency  domain  results  of  each  of  the  methods  described  above. 
For  the  time  domain,  we  plot  the  DRF  used  in  the  convolution  definition  of  the  polarization  from 
equation  3.  While  the  homogenization  method  used  here  directly  produces  a  dielectric  response 
function  defined  at  specified  time  points,  the  inverse  problem  approach  is  designed  to  simply  output 
a  single  parameter  set  for  use  in  a  Debye  polarization  model  formulation.  To  compare  the  two 
methods  in  the  time  domain  we  may  use  the  fact  that  the  Debye  model  DRF  is  given  by  (4). 
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In  order  to  examine  the  differences  in  the  frequency  domain  we  must  compare  the  resulting 
complex  permittivity.  For  the  Debye  polarization  model,  the  (relative)  complex  permittivity  is 
given  by 


e(w) 


eoo  + 


£s  Coo 
1  +  i(jjr 


However,  the  homogenization  approach  used  here  does  not  output  a  complex  permittivity  directly. 
One  may  take  the  DFT  of  the  DRF  to  get  the  permittivity.  This  presents  some  computation 
issues  when  attempting  to  plot  the  permittivity  for  the  entire  range  of  frequencies  of  interest. 
Therefore,  plots  involving  the  homogenization  method  are  incomplete  in  some  places,  and  care 
must  be  taken  in  interpreting  the  results.  See  Section  7  for  a  discussion  on  the  convergence  of  the 
permittivity  curves  as  the  discretization  parameters  are  reduced.  In  general,  the  homogenization 
results  tend  to  match  those  of  the  classical  Maxwell-Garnett  mixing  formula.  A  more  efficient 
approach  for  frequency  domain  comparison  would  be  to  compute  the  effective  permittivity  from 
the  homogenization  method  within  a  frequency  domain  formulation,  as  described  in  [6] .  However, 
the  main  purpose  of  this  effort  is  time  domain  comparisons,  therefore  we  only  present  the  frequency 
domain  results  to  emphasize  the  difference  between  the  methods.  Plots  of  the  real  and  imaginary 
parts  of  the  complex  permittivity  are  plotted  in  the  following  subsections. 


An  additional  depiction  of  the  complex  permittivity  is  the  Cole-Cole  diagram.  Cole  and  Cole 
[9]  proposed  an  empirical  function  for  use  in  fitting  dielectric  data.  The  Debye  model  presented 
in  Section  2.1  is  a  special  case  of  this  Cole-Cole  model.  The  Cole-Cole  function  (including  Debye) 
leads  to  a  distinctive  semicircular  plot  in  the  complex  permittivity  plane.  As  the  frequency  is 
varied,  a  plot  of  the  imaginary  part  of  the  complex  permittivity  versus  the  real  part  describes  a 
circular  arc  [11]  which  among  other  names  is  called  the  Cole-Cole  diagram. 


6.1  80%  —  20%  Volume  Fraction 

The  left  plot  of  Figure  9  depicts  the  DRF’s  for  the  inverse  problem  result,  the  solution  from  the 
homogenization  approach,  as  well  as  the  curves  for  a  Debye  model  of  pure  water  and  for  pure 
ethanol.  In  the  right  plot  of  Figure  9  we  compare  the  DRF  from  the  inverse  problem  and  homog¬ 
enization  approaches  to  the  MG  mixing  rule  result  and  the  DRF  obtained  by  using  the  weighted 
average  of  values  in  the  Debye  model  (4).  The  inverse  problem  DRF  is  nearly  indistinguishable 
from  that  of  the  weighted  averages;  however,  the  difference  between  the  corresponding  outputs 
from  the  interrogation  simulations  was  clearly  seen  in  Figure  6. 

The  real  and  imaginary  parts  of  the  permittivity  resulting  from  the  homogenization  approach 
and  the  inverse  problem  approach,  are  plotted  in  Figure  10,  along  with  the  values  for  pure  water  and 
pure  ethanol.  Both  methods  have  curves  that  lie  between  those  of  the  water  and  ethanol,  although 
more  closely  aligned  with  water.  However  the  main  difference  is  that  the  homogenization  approach 
has  a  peak  in  the  imaginary  part  of  the  permittivity  very  close  to  the  relaxation  frequency  of  water, 
whereas  the  inverse  problem  formulation  seems  to  have  a  peak  close  to  the  weighted  average  of 
relaxation  frequencies. 

In  Figure  11  the  real  and  imaginary  parts  of  the  permittivity  from  the  homogenization  approach 
and  the  inverse  problem  approach  are  compared  to  the  values  resulting  from  taking  the  weighted 
average  of  parameters  and  the  Maxwell-Garnett  mixing  method. 

The  peaks  in  the  plots  of  the  imaginary  parts  denote  the  frequency  at  which  the  maximum 
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Figure  9:  The  left  side  presents  a  plot  of  the  DRF  resulting  from  the  homogenization  approach  and 
the  inverse  problem  approach,  using  a  volume  fraction  of  80%  —  20%,  along  with  the  plots  for  the 
Debye  model  of  water  and  ethanol.  The  right  side  depicts  the  same  for  the  homogenization  and 
inverse  problem  approaches,  along  with  the  DRF’s  computed  using  the  Maxwell-Garnett  mixing 
formula  and  a  weighted  average.  The  DRF  of  the  weighted  average  approach  and  the  inverse 
problem  method  are  nearly  indistinguishable. 


Figure  10:  The  left  (right)  graph  depicts  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach  using  a  volume  fraction  of 
80%  —  20%,  along  with  the  values  for  pure  water  and  pure  ethanol. 


attenuation  of  the  field  in  the  material  occurs,  (in  the  Debye  case  it  is  known  as  the  relaxation 
frequency).  Note  that  both  the  Maxwell-Garnett  and  the  homogenization  methods  have  peaks  very 
close  to  that  of  water.  Both  the  parameter  estimation  method  and  the  approach  using  weighted 
average  of  parameters  result  in  peaks  very  close  to  the  weighted  average  of  the  relaxation  frequencies 
of  water  and  ethanol. 
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Figure  11:  The  left  (right)  graph  presents  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach,  along  with  the  values 
computed  using  a  weighted  average  and  the  M  ax  we  11-Garnett  mixing  formula,  using  a  volume 
fraction  of  80%  —  20%. 


Figure  12:  The  left  plot  depicts  the  Cole-Cole  diagram  of  the  complex  permittivity  resulting  from 
the  homogenization  approach  and  the  inverse  problem  approach,  along  with  the  values  from  the 
Debye  model  for  water  and  ethanol  using  a  volume  fraction  of  80%  —  20%.  The  right  plot  presents 
the  same  for  the  homogenization  and  inverse  problem  approaches,  along  with  the  values  computed 
using  the  Maxwell-Garnett  mixing  formula  and  a  weighted  average. 


The  Cole-Cole  diagrams  corresponding  to  the  complex  permittivities  shown  in  Figures  10  and 
11  are  shown  in  the  left  and  right  plots  of  Figure  12. 

The  Cole-Cole  diagram  for  a  Debye  model  is  a  semi-circle  and  the  peaks  correspond  to  the 
relaxation  frequency.  As  both  the  weighted  average  values  and  the  inverse  problem  approach 
use  the  Debye  model,  their  diagrams  are  semi-circles.  Both  the  Maxwell-Garnett  rule  and  the 
homogenization  method  have  more  complex  Cole-Cole  diagrams  as  evidenced  by  the  bend  in  the 
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curves.  This  discrepancy  from  the  Debye  model  is  more  pronounced  at  the  lower  frequencies  for 
this  volume  fraction. 


6.2  50%  —  50%  Volume  Fraction 

The  left  plot  of  Figure  13  depicts  the  DRF’s  for  the  inverse  problem  result,  the  solution  from 
the  homogenization  approach,  as  well  as  the  curves  for  a  Debye  model  of  pure  water  and  for 
pure  ethanol.  In  the  right  plot  of  Figure  13  we  compare  the  DRF  from  the  inverse  problem  and 
homogenization  approaches  to  the  MG  mixing  rule  result  and  the  DRF  obtained  by  using  the 
weighted  average  of  values  in  the  Debye  model.  The  inverse  problem  DRF  is  very  similar  to  that 
of  the  weighted  averages;  however,  the  difference  between  the  corresponding  outputs  from  the 
interrogation  simulations  was  clearly  seen  in  Figure  7. 


Figure  13:  The  left  side  presents  a  plot  of  the  DRF  resulting  from  the  homogenization  approach 
and  the  inverse  problem  approach,  along  with  the  plots  for  the  Debye  model  of  water  and  ethanol, 
using  a  volume  fraction  of  50%  —  50%.  The  right  side  depicts  the  same  for  the  homogenization  and 
inverse  problem  approaches,  along  with  the  DRF’s  computed  using  the  Maxwell-Garnett  mixing 
formula  and  a  weighted  average.  The  DRF  of  the  weighted  average  approach  and  the  inverse 
problem  method  are  nearly  indistinguishable. 


The  real  and  imaginary  parts  of  the  permittivity  resulting  from  the  homogenization  approach 
and  the  inverse  problem  approach,  are  plotted  in  Figure  14,  along  with  the  values  for  pure  water  and 
pure  ethanol.  In  Figure  15  the  real  and  imaginary  parts  of  the  permittivity  from  the  homogenization 
approach  and  the  inverse  problem  approach  are  compared  to  the  values  resulting  from  taking  the 
weighted  average  of  parameters  and  the  Maxwell-Garnett  mixing  method. 

The  peaks  in  the  plots  of  the  imaginary  parts  denote  the  frequency  at  which  the  maximum 
attenuation  of  the  field  in  the  material  occurs,  (in  the  Debye  case  it  is  known  as  the  relaxation 
frequency).  Note  that  both  the  Maxwell-Garnett  and  the  homogenization  methods  still  have  peaks 
relatively  close  to  that  of  water,  as  compared  to  ethanol,  considering  that  it  is  a  50%  —  50%  mixture. 
Both  the  parameter  estimation  method  and  the  approach  using  weighted  average  of  parameters 
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Figure  14:  The  left  (right)  graph  depicts  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach,  using  a  volume  fraction 
of  50%  —  50%,  along  with  the  values  for  pure  water  and  pure  ethanol. 


Figure  15:  The  left  (right)  graph  presents  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach,  along  with  the  values 
computed  using  a  weighted  average  and  the  Maxwell-Garnett  mixing  formula,  using  a  volume 
fraction  of  50%  —  50%. 


result  in  peaks  very  close  to  the  average  of  the  relaxation  frequencies  of  water  and  ethanol. 

The  Cole-Cole  diagrams  corresponding  to  the  complex  permittivities  shown  in  Figures  14  and 
15  are  shown  in  the  left  and  right  plots  of  Figure  16.  As  both  the  weighted  average  values  and 
the  inverse  problem  approach  use  the  Debye  model,  their  diagrams  are  again  semi-circles.  The 
Maxwell-Garnett  rule  and  the  homogenization  method  are  again  more  complicated,  and  have  an 
even  more  pronounced  bend  in  the  curve  than  the  80%  —  20%  case.  The  frequency  at  which  is 
occurs  is  also  closer  to  the  peak  frequency  than  the  previous  case. 
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Figure  16:  The  left  plot  depicts  the  Cole-Cole  diagram  of  the  complex  permittivity  resulting  from 
the  homogenization  approach  and  the  inverse  problem  approach,  along  with  the  values  from  the 
Debye  model  for  water  and  ethanol,  using  a  volume  fraction  of  50%  —  50%.  The  right  plot  presents 
the  same  for  the  homogenization  and  inverse  problem  approaches,  along  with  the  values  computed 
using  the  Maxwell-Garnett  mixing  formula  and  a  weighted  average. 


6.3  20%  —  80%  Volume  Fraction 

The  left  plot  of  Figure  17  depicts  the  DRF’s  for  the  inverse  problem  result,  the  solution  from 
the  homogenization  approach,  as  well  as  the  curves  for  a  Debye  model  of  pure  water  and  for 
pure  ethanol.  In  the  right  plot  of  Figure  17  we  compare  the  DRF  from  the  inverse  problem  and 
homogenization  approaches  to  the  MG  mixing  rule  result  and  the  DRF  obtained  by  using  the 
weighted  average  of  values  in  the  Debye  model.  The  inverse  problem  DRF  is  very  similar  to  that 
of  the  weighted  averages;  however,  the  difference  between  the  corresponding  outputs  from  the 
interrogation  simulations  was  clearly  seen  in  Figure  8. 

The  real  and  imaginary  parts  of  the  permittivity  resulting  from  the  homogenization  approach 
and  the  inverse  problem  approach,  are  plotted  in  Figure  18,  along  with  the  values  for  pure  water  and 
pure  ethanol.  In  Figure  19  the  real  and  imaginary  parts  of  the  permittivity  from  the  homogenization 
approach  and  the  inverse  problem  approach  are  compared  to  the  values  resulting  from  taking  the 
weighted  average  of  parameters  and  the  Maxwell-Garnett  mixing  method. 

The  peaks  in  the  plots  of  the  imaginary  parts  denote  the  frequency  at  which  the  maximum 
attenuation  of  the  field  in  the  material  occurs,  (in  the  Debye  case  it  is  known  as  the  relaxation 
frequency).  Note  that  both  the  Maxwell-Garnett  and  the  homogenization  methods  still  have  peaks 
relatively  close  to  that  of  water,  as  compared  to  ethanol,  considering  that  it  is  a  20%  —  80%  mixture. 
Both  the  parameter  estimation  method  and  the  approach  using  weighted  average  of  parameters 
result  in  peaks  very  close  to  the  average  of  the  relaxation  frequencies  of  water  and  ethanol. 

The  Cole-Cole  diagrams  corresponding  to  the  complex  permittivities  shown  in  Figures  18  and 
19  are  shown  in  the  left  and  right  plots  of  Figure  20.  As  both  the  weighted  average  values  and 
the  inverse  problem  approach  use  the  Debye  model,  their  diagrams  are  again  semi-circles.  The 
Maxwell-Garnett  rule  and  the  homogenization  method  are  again  more  complicated,  and  have  an 
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Figure  17:  The  left  side  presents  a  plot  of  the  DRF  resulting  from  the  homogenization  approach 
and  the  inverse  problem  approach,  along  with  the  plots  for  the  Debye  model  of  water  and  ethanol, 
using  a  volume  fraction  of  20%  —  80%.  The  right  side  depicts  the  same  for  the  homogenization  and 
inverse  problem  approaches,  along  with  the  DRF’s  computed  using  the  M axwell- G ar net t  mixing 
formula  and  a  weighted  average.  The  DRF  of  the  weighted  average  approach  and  the  inverse 
problem  method  are  nearly  indistinguishable. 


Figure  18:  The  left  (right)  graph  depicts  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach,  using  a  volume  fraction 
of  20%  —  80%,  along  with  the  values  for  pure  water  and  pure  ethanol. 


even  more  pronounced  bend  in  the  curve  than  the  80%  —  20%  case.  The  frequency  at  which  is 
occurs  is  also  closer  to  the  peak  frequency  than  the  previous  case. 
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Figure  19:  The  left  (right)  graph  presents  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  homogenized  DRF  and  the  inverse  problem  approach,  along  with  the  values 
computed  using  a  weighted  average  and  the  M  ax  we  11-Garnett  mixing  formula,  using  a  volume 
fraction  of  20%  —  80%. 


Figure  20:  The  left  plot  depicts  the  Cole-Cole  diagram  of  the  complex  permittivity  resulting  from 
the  homogenization  approach  and  the  inverse  problem  approach,  along  with  the  values  from  the 
Debye  model  for  water  and  ethanol,  using  a  volume  fraction  of  20%  —  80%.  The  right  plot  presents 
the  same  for  the  homogenization  and  inverse  problem  approaches,  along  with  the  values  computed 
using  the  Maxwell-Garnett  mixing  formula  and  a  weighted  average. 


7  Convergence  of  DFT  for  Homogenization  DRF 


In  order  to  display  the  complex  permittivity  resulting  from  the  DRF  computed  with  the  homoge¬ 
nization  approach,  one  must  take  the  DFT  of  the  DRF  at  many  time  points.  Further,  in  order  for 
the  maximum  frequency  component  to  be  large  enough  to  display  useful  information,  the  time  step 
must  be  very  small.  For  example,  in  order  to  compute  the  complex  permittivity  up  to  /  =  5  x  1012 
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Hz  (where  the  real  part  of  the  permittivity  converges  to  Coo),  the  time  step  for  calculating  the  DRF 
in  the  homogenization  method  must  be  10“13  seconds.  Further,  for  sufficient  resolution  down  to 
/  =  108  Hz  (where  the  real  part  of  the  permittivity  converges  to  es),  the  number  of  time  steps 
must  be  greater  than  20000.  In  our  numerical  experiments,  each  time  step  takes  on  the  order  of  .1 
seconds,  as  a  linear  system  for  the  finite  element  method  must  be  solved.  The  following  subsections 
show  computations  of  the  DFT  for  various  DRF’s  computed  with  successively  smaller  time  steps. 


7.1  80%  —  20%  Volume  Fraction 

We  compute  the  homogenized  DRF  for  the  volume  fraction  80%  —  20%  water-ethanol  using  time 
steps  dt  =  { 10  14, 10~13, 10~12}.  The  resulting  curves  are  shown  in  the  left  plot  of  Figure  21.  The 
convergence  in  the  time  domain  is  quickly  apparent,  however,  many  more  steps  would  be  required 
for  full  convergence  of  the  frequency  domain  results,  as  seen  in  the  right  plot  which  contains  the 
Cole-Cole  diagram  for  each  time  step.  The  real  and  imaginary  parts  of  the  permittivity  resulting 


Figure  21:  The  left  side  presents  a  plot  of  the  DRF  resulting  from  the  homogenization  approach 
using  time  steps  dt  =  {10-14, 10  13, 10-12}.  for  the  case  of  an  80%  —  20%  volume  fraction  of  water 
and  ethanol.  The  right  plot  depicts  the  Cole-Cole  diagram  of  the  complex  permittivity  resulting 
from  the  homogenization  approach  for  each  time  step.  The  incompleteness  is  due  to  lack  of  sufficient 
time  steps.  105  steps  are  shown. 


from  the  homogenization  approach  using  time  steps  dt  =  {10~14, 10~13, 10-12}  are  plotted  in  Figure 
22.  The  high  frequency  part  of  the  curves  are  complete  and  appear  to  be  converging  to  a  limiting 
effective  complex  permittivity. 


7.2  20%  —  80%  Volume  Fraction 

For  the  case  of  a  volume  fraction  of  20%  —  80%  water-ethanol  we  use  the  time  steps  dt  =  {10~16, 
10_1  ,  10  13}.  The  resulting  curves  are  shown  in  the  left  plot  of  Figure  23.  Although  fewer  time 
steps  were  computed  in  the  dt  =  10“16  case,  clear  agreement  in  the  time  domain  is  reached  with 
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Figure  22:  The  left  (right)  graph  depicts  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  DFT  of  the  homogenized  DRF,  using  a  volume  fraction  of  80%  —  20%,  computed 
with  time  steps  dt  =  {10~14, 10~13, 10~12}. 


only  dt  =  10~14.  The  right  plot  of  Figure  23  depicts  the  Cole-Cole  diagram  for  each  time  step.  All 
three  resolutions  agree  at  the  higher  frequencies  (the  left  side  of  the  plot),  which  is  where  the  effect 
of  the  composite  material  is  most  obvious.  The  overlapping  circles  is  characteristic  of  mixtures  of 
Debye  substances,  and  the  diagrams  here  agree  nicely  with  the  classical  Maxwell-Garnett  mixing 
rule  shown  in  Figure  20.  The  real  and  imaginary  parts  of  the  permittivity  resulting  from  the 


Figure  23:  The  left  side  presents  a  plot  of  the  DRF  resulting  from  the  homogenization  approach 
using  time  steps  dt  =  {10-16, 10~14, 10  13}.  for  the  case  of  an  20%  —  80%  volume  fraction  of  water 
and  ethanol.  The  right  plot  depicts  the  Cole-Cole  diagram  of  the  complex  permittivity  resulting 
from  the  homogenization  approach  for  each  time  step.  The  inaccuracies  due  to  the  incompleteness 
of  the  time  domain  DRF  are  clearly  apparent. 

homogenization  approach  using  time  steps  dt  =  {10~16, 10-14, 10~13},  are  plotted  in  Figure  24. 
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The  high  frequency  part  of  the  curves  appear  to  coincide  better  than  for  the  80%  —  20%  case, 
however  the  lower  frequencies  are  not  sufficiently  resolved  to  be  able  to  estimate  the  effective  static 
permittivity  es. 
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Figure  24:  The  left  (right)  graph  depicts  the  real  (negative  imaginary)  part  of  the  permittivity 
resulting  from  the  DFT  of  the  homogenized  DRF,  using  a  volume  fraction  of  20%  —  80%,  computed 
with  time  steps  dt  =  {10-16, 1CT14, 1CT13}. 


8  Conclusions 


The  Maxwell-Garnett  mixing  rule  is  only  valid  for  spherical  or  circular  inclusions  (ellipsoidal  inclu¬ 
sions  can  also  be  considered).  The  periodic  unfolding  homogenization  method  allows  for  any  two  or 
three  dimensional  geometry.  It  is  important  to  note,  however,  that  the  results  of  all  homogenization 
techniques  (including  MG)  are  valid  when  the  size  of  the  inclusion  are  small  in  comparison  with 
the  wave  length  of  the  sources  and  fields.  For  higher  frequencies  scattering  effects  are  no  longer 
negligible.  A  second  advantage  of  the  homogenization  method  presented  here  is  the  flexibility  it 
affords  in  assumptions  about  material  polarization  laws.  In  this  paper  we  have  computed  the  DRF 
for  a  Debye-Debye  mixture  using  a  recursive  convolution  method.  The  recursive  convolution  ap¬ 
proach  can  be  extended  to  the  cases  of  mixtures  involving  Lorentz  or  higher  order  dispersive  media 
in  which  the  DRF  is  exponential  in  nature. 

A  mixture  of  materials,  which  may  each  be  described  by  the  same  simple  model,  in  general 
will  not  be  accurately  described  by  that  same  model.  For  example,  a  Debye-Debye  mixture  has 
an  effective  complex  permittivity  which  is  not  of  Debye  form,  nor  even  of  Lorentz  form.  The 
homogenization  methods  described  here  attempt  to  capture  this  complexity.  However,  for  certain 
applications  it  is  not  necessary  to  have  a  high  level  of  accuracy  for  those  frequencies  which  do  not 
appear  in  the  problem,  especially  if  the  added  complexity  necessary  for  increased  accuracy  is  cost- 
prohibitive.  For  example,  in  a  simulation  with  an  interrogating  signal  of  narrow  frequency  band, 
the  complex  permittivity  need  only  be  accurate  in  those  frequencies.  Using  an  inverse  problem 
formulation  designed  around  the  specific  fields  of  interest  can  give  a  cost-effective  method  for 
modeling  complex  materials  with  relatively  simple  formulas. 
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As  mentioned  earlier,  for  both  the  homogenization  method  based  on  periodic  unfolding  as  well 
as  the  Maxwell-Garnett  model,  knowledge  of  the  volume  proportions  of  the  components  of  the 
mixture  is  needed  in  order  to  calculate  effective  parameters  of  this  composite  dielectric.  However, 
the  parameter  estimation  method  is  capable  of  using  experimental  data  in  the  inverse  problem,  in 
which  case  there  would  be  no  assumptions  needed  on  the  volume  fraction  or  dielectric  parameters 
except  perhaps  in  formulating  a  sufficiently  accurate  initial  value  for  faster  convergence  of  the 
optimization  routine.  In  generating  synthetic  data  however,  one  must  of  course  have  the  volume 
proportions,  as  well  as  dielectric  parameters  for  each  material,  including  possibly  distributions  of 
parameters.  While  there  is  a  large  initial  investment  in  computational  time  to  perform  the  many 
forward  solves  needed  to  determine  the  optimal  parameter  set,  having  effective  parameters  for  a 
simplified  model  which  closely  matches  the  data  of  interest  may  lead  to  significant  savings  in  later 
simulations  such  as  in  geometric  inverse  problems  involving  crack  or  flaw  detection.  One  of  the 
biggest  advantages  to  using  the  inverse  problem  parameter  estimation  approach  is  that  it  may  be 
tailored  to  specific  applications  with  narrow  frequency  bands  such  as  those  that  often  occur  in  the 
fields  of  non-destructive  evaluation  (NDE)  or  remote  sensing. 
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